!
! Conversion between astronomical coordinates
!
! $Id$
!
! --------------------------------------------------------------------

Module traco

  implicit none

  integer, parameter, private :: db = selected_real_kind(15)
  real(db), parameter, private :: RAD =  57.29577951308232286464772_db

contains

  function gmst(jd)

    implicit none
    real(db) :: gmst
    real(db), intent(in) :: jd
    

!  Greenwich siderical time
!
!  jd is a full Julian date
!
!  The precision is better than 1 second. 
!  According to Astronomical Almanac 2000.


    real(db) :: tu,t

    tu = (jd - 2451545.0_db)/36525.0_db
    t = 24110.54841_db + tu*(8640184.812866_db + tu*(0.093104_db-6.2e-6_db*tu))
    gmst = mod(t/3600.0_db + 24.0_db*(jd - aint(jd)) + 12.0_db,24.0_db)

  end function gmst


  function lmst(jd,longitude)

    implicit none
    real(db) :: lmst
    real(db), intent(in) :: jd, longitude


!  local siderical time
!
!  jd is a full Julian date
!  lamda is a longitude in degrees, -west... +east
!
!  The precision is better than 1 second. 
!  According to Astronomical Almanac 2000.


    lmst = mod(gmst(jd) + longitude/15.0_db,24.0_db)

  end function lmst



  function hangle(lmst,ra)

    implicit none
    real(db) :: hangle
    real(db), intent(in) :: lmst, ra


!   hour angle


    hangle = mod(lmst - ra,360.0_db)

  end function hangle

  
  subroutine eq2hor(ha, dec, latitude, az, elev)

    implicit none
    real(db), intent(in) :: ha,dec,latitude
    real(db), intent(out) :: az, elev


!   
! equatorial to horizontal coordinates
!
!  all arguments in degrees
!

    real(db) :: sinh, cosh, sind, cosd, sinl, cosl, x,y,z,r
    
    sinh = sin(ha/RAD)
    cosh = cos(ha/RAD)
    sind = sin(dec/RAD)
    cosd = cos(dec/RAD)
    sinl = sin(latitude/RAD)
    cosl = cos(latitude/RAD)

    x = -cosh*cosd*sinl + sind*cosl
    y = -sinh*cosd
    z = cosh*cosd*cosl + sind*sinl

    r = sqrt(x**2 + y**2)
    if( abs(r) > epsilon(r) )then
       az = RAD*atan2(y,x)
    else
       az = 0_db
    end if
    if( az < 0_db ) az = az + 360.0_db
    elev = RAD*atan2(z,r)

  end subroutine eq2hor


  subroutine hor2eq(az, elev, latitude, ha,dec)

    implicit none
    real(db), intent(in) :: az,elev,latitude
    real(db), intent(out) :: ha, dec

!
!  horizontal to equatorial coordinates
!
!  all arguments in degrees
!
    
    real(db) :: sina, cosa, sine, cose, sinl, cosl, x, y, z, r

    sina = sin(az/RAD)
    cosa = cos(az/RAD)
    sine = sin(elev/RAD)
    cose = cos(elev/RAD)
    sinl = sin(latitude/RAD)
    cosl = cos(latitude/RAD)

    x = -cosa*cose*sinl + sine*cosl
    y = -sina*cose
    z = cosa*cose*cosl + sine*sinl
    
    r = sqrt(x**2 + y**2)
    if( abs(r) > epsilon(r) )then
       ha = RAD*atan2(y,x)
    else
       ha = 0.0
    endif
    dec = RAD*atan2(z,r)

  end subroutine hor2eq

  
  function refract(z)

    implicit none
    real(db) :: refract
    real(db), intent(in) :: z

!
!     compute refraction angle in degrees 
!
!    Smart: Textbook on spherical astronomy
!
!    constants for pressure 760mmHg, 10deg C with 
!    suffucient accuracy for z < 75 deg 
!

    real(db) :: tanz

    tanz = tan(z/RAD)
    refract = (58.16_db*tanz - 0.067_db*tanz*tanz*tanz)/3600.0_db

  end function refract


  function airmass(z)

    implicit none
    real(db) :: airmass
    real(db), intent(in) :: z

!
!     compute airmass, 
!
!      young&irvine: aj,72,945,(1967)

    real(db) :: secz

    secz = 1.0/cos(z/RAD)
    airmass = secz*(1.0_db - 1.2e-3_db*(secz**2 - 1.0_db))

  end function airmass

end Module traco

